ADIM: Bayesian GLM. MCMC methods

.title[
# ADIM: Bayesian GLM. MCMC methods
]
.subtitle[
## Master’s Degree in Data Analysis, Process Improvement and Decision Support Engineering
]
.author[
### </br></br> <font size="6"> Joaquín Martínez-Minaya, 2024-12-02 </font> </br> </br> <a href="http://vabar.es/" style="color:white;"> VAlencia BAyesian Research Group </a> </br> <a href="https://smeg-bayes.org/ " style="color:white;"> Statistical Modeling Ecology Group </a> </br> <a href="https://giem.blogs.upv.es/" style="color:white;"> Grupo de Ingeniería Estadística Multivariante </a> </br>
]
.date[
### <a href="jmarmin@eio.upv.es" style="color:white;"> jmarmin@eio.upv.es </a>
]

---

# Motivation example
---
# Heart Disease

- the .hlb[myocardial infarction] (MI): `$y = 1$` if MI occurrence, or `$y = 0$` if No MI occurrence; and
  - .hlb[Age60]: Patients aged `$\geq 60$` (1) versus `$< 60$` (0).
  - .hlb[Systolic blood pressure (SBP140)]: SBP `$\geq 140$` mmHg (1) versus `$< 140$` mmHg (0).

- .hlbred[Objective]:
  - Evaluate the association of age60 and sbp140 with MI probability.
  - Interpret the odds ratio (OR) for both predictors.
]

Table: Summary of Data from Study

| y | age60 | sbp140 |
|:-:|:-----:|:------:|
| 0 |  <60  | >=140  |
| 0 | >=60  |  <140  |
| 0 |  <60  | >=140  |
| 0 | >=60  | >=140  |
| 0 | >=60  |  <140  |
| 1 |  <60  |  <140  |
]
---

# Bayesian Logistic Regression Model

- .hlbpurple[Logistic regression] is used to model MI's probability based on age60 and sbp140.

- .hlb[Likelihood]
  $$ y_i  \sim  \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$
  using logit link: 
  
    $$ logit(\pi_i)  =  \beta_0 + \beta_1 age60_i + \beta_2 sbp140_i$$ 
  
- .hlb[Prior distributions] (weakly-informative):
    \begin{eqnarray}
      \beta_0 & \sim & \mathcal{N}(0, 10^3), \nonumber \\
      \beta_1 & \sim & \mathcal{N}(0, 10^3), \nonumber \\
      \beta_2 & \sim & \mathcal{N}(0, 10^3), \nonumber
    \end{eqnarray}
    
    
--
.hlbred[Note: There are no conjugate priors available for the logistic regression model.]

---

#Table of contents

## 1. Bayesian computation. MCMC methods

## 2. Bayesian Software for MCMC

---
class: center, middle, animated, rotateInUpRight, inverse
# 1. Bayesian computation. MCMC methods

---

# Monte Carlo Methods

- Draw .hlb[realizations of a random variable] for which only its density function is (fully or partially) known.

``` r
  x <- rnorm(1000, mean = 1, sd = 2)
```

![](data:image/png;base64,#S1_2_bayesian_computation_files/figure-html/unnamed-chunk-9-1.png)
]

.hlbpurple[Monte Carlo Integration]
- Computing the mean of a N(1, 2), 
  - `$E(X) = \int_{-\infty}^{\infty} x \cdot f(x) dx$`
- Using .hlb[Monte Carlo integration]:
  - Simulate from N(1, 2^2): `$\phi^1, \ldots, \phi^N$`.
  - Compute the mean of the simulated values: `$E(X) \approx \frac{1}{N} \sum_{i=1}^N \phi^i$`
- Doing `summary` of the simulation, we compute more measures:

|   Min.| 1st Qu.| Median| Mean| 3rd Qu.|  Max.|
|------:|-------:|------:|----:|-------:|-----:|
| -6.065|  -0.334|      1|    1|   2.332| 8.086|
]

---

# Markov Chain Monte Carlo

- A Markov chain is a .hlbred[stochastic sequence of numbers] where each value in the sequence depends only upon the last.

- If `$\phi^1, \phi^2, \ldots, \phi^N$` is a sequence of numbers, then `$\phi^2$` is only a function of `$\phi^1$`, `$\phi^3$` of `$\phi^2$`, etc.

- Under certain conditions, the distribution over the states of the .hlb[Markov chain] will .hlb[converge to a stationary distribution].

- The .hlb[stationary distribution is independent of the initial starting values] specified for the chains.

- AIM: construct a Markov chain such that .hlb[the stationary distribution is equal to the posterior distribution] `$p(\theta \mid x)$`.

- We combine Markov Chain with Monte Carlo simulation --> .hlbpurple[Markov chain Monte Carlo (MCMC)].

- They were proposed by first time in the Statistics area by <a href = "https://www.tandfonline.com/doi/abs/10.1080/01621459.1990.10476213"> Gelfand and Smith (1990) </a>.

---
# Posterior distribution
## Estimating the probability to score a penalty

### - .hlb[Likelihood]
$$p(\boldsymbol{y} \mid \pi) =  \pi^{k}(1-\pi)^{N-k} $$

### - .hlb[Prior distribution]
$$p(\pi) = \pi^{a-1}(1-\pi)^{b-1} $$

### - .hlb[Posterior distribution]
`$$p(\pi \mid \boldsymbol{y})  \propto  p(\boldsymbol{y} \mid \pi) \times p(\pi) \propto \pi^{k + a - 1}(1 - \pi)^{N - k + b - 1} = p^*(\pi)$$`
---
# MCMC: Metropolis-Hastings (MH)

1. Starting value `$\pi^{(0)}$`

2. For `$t = 1, \ldots, T$`
  - .hlb[We define a proposal distribution] (Usually similar to the objective distribution). In this case, `$q(\pi \mid \pi^{(t-1)}) \sim logit-N(\pi^{(t-1)}, \sigma = 0.5)$`. .hlb[Simulate] `$\pi^{(prop)}$` from it.
  - Compute .hlbred[probability of acceptance]:
    `$$\alpha = \text{min}\left(1, \frac{p^*(\pi^{(prop)})q(\pi^{(t-1)} \mid \pi^{(prop)})}{p^*(\pi^{(t-1)})q(\pi^{(prop)} \mid \pi^{(t-1)})}\right)$$`
  - Generate a .hlb[random number] `$u$` from the Uniform(0, 1).
     - `$\pi^{(t+1)} = \pi^{(prop)}$`, if `$u \geq \alpha$`, 
     - `$\pi^{(t+1)} = \pi^{(t)}$`, if `$u< \alpha$`
      
3. Finally, we .hlb[obtain] `$\pi^{0}, \pi^{1}, \ldots, \pi^{T}$` which is .hlb[a simulation of the posterior distribution].

---
# MCMC: Metropolis-Hastings (MH)
.center[
![](data:image/png;base64,#images/mcmc.png)
]

---
# Approaching probability of score using MH

- Play the video from .hlb[minute 4:44].

]

![](data:image/png;base64,#images/sim_acf_chain1.png)
]

---
# MCMC. Burnin and thin
.left-column9[
![](data:image/png;base64,#images/sims_chain_pi.png)
]

.right-column9[
![](data:image/png;base64,#images/acf_sim2.png)
![](data:image/png;base64,#images/posterior_pi_sim.png)

]

---

# 2. Bayesian Software for MCMC

---

# Software

- <a href="https://mc-stan.org/"> Stan </a>

- <a href = "https://paulbuerkner.com/brms/"> brms </a>: allow for easy Bayesian Inference using Hamiltonian Monte Carlo
- <a href="https://www.r-inla.org/"> INLA </a> (it does not use MCMC methods, but it is a very powerful tool)

- <a href = "https://inlabru-org.github.io/inlabru/"> inlabru </a>: facilitate Bayesian Inference in Spatio-temporal models
  
- <a href= "https://cran.r-project.org/web/packages/MCMCpack/index.html"> MCMC pack </a> (With this R-package, you can use MCMC methods with similar notation as usually use in R)
  
- <a href="https://r-nimble.org/"> Nimble </a>

]

---

# Bayesian Logistic Regression using JAGS
.left-column8[
- .hlb[Likelihood]
  $$ y_i  \sim  \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$
  using logit link: 
  
    $$ logit(\pi_i)  =  \beta_0 + \beta_1 age60_i + \beta_2 sbp140_i$$ 
  
- .hlb[Prior distributions] (weakly-informative):
      `$$\beta_0  \sim  \mathcal{N}(0, 10^3),$$` 
      `$$\beta_1  \sim  \mathcal{N}(0, 10^3),$$` 
      `$$\beta_2  \sim  \mathcal{N}(0, 10^3).$$` 
]
.right-column9[

``` r
model_string <- "
model {
  for (i in 1:N) {
    y[i] ~ dbern(pi[i])
    logit(pi[i]) <- beta0 + 
      beta1 * age60[i] + 
      beta2 * sbp140[i]
  }
  # Priors for regression coefficients
  beta0 ~ dnorm(0, 0.001)
  beta1 ~ dnorm(0, 0.001)  
  beta2 ~ dnorm(0, 0.001)  
}"
```

Check `S1-JAGS-heart_attack.Rmd` for the complete solution
]
---

# Bayesian Logistic Regression using brms

.left-column8[
- .hlb[Likelihood]
  $$ y_i  \sim  \text{Bernoulli}(\pi_i)\,, i= 1, \ldots, 400\,, $$
  using logit link: 
  
  $$ \text{logit}(\pi_i)  =  \beta_0 + \beta_1 \, \text{age60}_i + \beta_2 \, \text{sbp140}_i $$ 
  
- .hlb[Prior distributions] (weakly-informative):
  `$$\beta_0  \sim  \mathcal{N}(0, 10^3),$$` 
  `$$\beta_1  \sim  \mathcal{N}(0, 10^3),$$` 
  `$$\beta_2  \sim  \mathcal{N}(0, 10^3).$$` 
]

``` r
formula <- bf(y ~ age60 + sbp140, family = bernoulli())

# Fit the model using brms
fit_brms <- brm(formula,
    data = data_hattack,
    prior = priors,
    chains = 3,      # Number of MCMC chains
    iter = 5000,     # Total number of iterations per chain
    warmup = 1000,   # Number of iterations for warm-up
    thin = 1,        # Thinning interval
    seed = 123,      # Seed for reproducibility
  )

# Summary of the fitted model
summary(fit)
```

Check `S1-brms-heart_attack.Rmd` for the complete solution
]

---
# What we have learned so far

- ### The biggest challenge for Bayesian inference has always been the .hlb[computational power] that more complex models require.

- ### .hlb[MCMC methods allow computationally estimating posterior distributions] that are analytically intractable.

- ### Today, there are many, many teams of researchers developing computational techniques for computing a posteriori distributions.

---

# References

- .hlb[Stan & R]: Lambert, B. (2018). A Student’s Guide to Bayesian Statistics. SAGE Publications.

- .hlb[JAGS]: Plummer, M. (2019). JAGS User Manual Version 4.3.0.

- .hlb[OpenBUGS]: Cowles, M. K. (2013). Applied Bayesian statistics: with R and OpenBUGS examples (Vol. 98). Springer Science & Business Media.

- .hlb[WinBUGS]: Ntzoufras, I. (2011). Bayesian modeling using WinBUGS. John Wiley & Sons.

- .hlb[Stan]: Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, Donald B. Rubin (2013). Bayesian Data Analysis. Chapman and Hall/CRC

- .hlb[INLA]: Gómez-Rubio, V. (2020). Bayesian inference with INLA. CRC Press.

---
#References

- http://wlm.userweb.mwn.de/R/wlmRcoda.htm

- https://rpubs.com/FJRubio/IntroMCMC

- https://darrenjw.wordpress.com/tag/mcmc/

- https://www.tweag.io/blog/2019-10-25-mcmc-intro1/

---
class: inverse, left, middle, animated, rotateInUpLeft
</br>
# ADIM: Bayesian GLM. MCMC methods
## Master's Degree in Data Analysis, Process Improvement and Decision Support Engineering

</br> <font size="6"> Joaquín Martínez-Minaya, 2024-10-21 </font> 
</br> 
<a href="http://vabar.es/" style="color:white;" > VAlencia BAyesian Research Group </a> </br> <a href="https://smeg-bayes.org/ " style="color:white;"> Statistical Modeling Ecology Group </a> </br> <a href="https://giem.blogs.upv.es/" style="color:white;"> Grupo de Ingeniería Estadística Multivariante </a> </br>

<a href="jmarmin@eio.upv.es" style="color:white;"> jmarmin@eio.upv.es </a>